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ABSTRACT 

We study the prospects for extracting detailed statistical properties of the neutral Hydrogen distri- 
bution during the era of reionization using the brightness temperature fluctuations from redshifted 21 
cm line emission. Detection of this signal is complicated by contamination from foreground sources such 
as diffuse Galactic synchrotron and free-free emission at low radio frequencies, extragalactic free-free 
emission from ionized regions and radio point sources. We model these foregrounds to determine the 
extent to which 21 cm fluctuations can be detected with upcoming experiments. We find that not only 
the level of correlation from one frequency to another, but also the functional form of the foreground 
correlations has a substantial impact on foreground removal. We calculate how well the angular power 
spectra of the 21cm fluctuations can be determined. We also show that the large-scale bias of the neutral 
Hydrogen gas distribution with respect to the density fleld, can be determined with high precision and 
used to distinguish between different reionization histories. 

Subject headings: cosmology: theory — large scale structure of Universe — diffuse radiation 



1. INTRODUCTION 

Recently there has been a great deal of interest in the 
reionization of the intergalactic medium (Cen 2003; Chen 
et al. 2003; Cen & McDonald 2002; Gnedin 2004; Wyithe 
& Loeb 2003; Haiman & Holder 2003; Sokasian et al. 2003; 
Miralda-Escude et al. 2000). This interest was sparked by 
the detection (Bennett et al. 2003; Kogut et al. 2003) of a 
correlation between polarization and temperature on large 
angular scales (Zaldarriaga 1997). The amplitude of the 
signal suggests that reionization began much earlier than 
z ~ 6, which is when the inter-galactic medium became 
completely ionized as we know from the detection of Gunn- 
Peterson troughs (Gunn & Peterson 1965; Becker et al. 
2001; Fan et al. 2002). 

Now that the ionization history appears to be quite rich, 
we wish to study the epoch in detail. Unfortunately, the 
CMB is mostly sensitive only to the integral history of 
reionization. Although different reionization histories with 
the same optical depth do give different large angular scale 
polarization power spectra, one cannot easily use these 
changes to fully reconstruct the reionization history as a 
function of redshift due to the large cosmic variance on 
large angular scales (Kaplinghat et al. 2003; Hu & Holder 
2003; Haiman & Holder 2003). 

Another way to study the reionization history of the 
Universe in more detail is through the observation of the 
distribution of neutral Hydrogen based on the spin-flip 
transition at a rest wavelength of 21 cm. Plans for up- 
coming low-frequency radio experiments such as PAST"'^ 

-"^http:/ /web. phys.cmu.edu/-past/ 



(Pen et al. 2004), Mileura Widefield Array (MWA)^, Low 
Frequency Array (LOFAR)'^ and Square Kilometer Array 
(SKA)"* have now motivated detailed study of the 21 cm 
background associated with neutral hydrogen during, and 
prior to, reionization (e.g., Scott & Rees 1990; Madau et al. 
1997; Gnedin & Ostriker 1997; Zaldarriaga et al. 2004; 
Morales & Hewitt 2003; Kumar et al. 1995; Tozzi et al. 
2000; Ihev et al. 2002; Bharadwaj & Ah 2005). Unlike the 
CMB, the advantage here is that by a priori selecting the 
observed frequency, one can probe the neutral content of 
the Universe at a given redshift, thereby obtaining a tomo- 
graphic view of the reionization process. This way, one can 
reconstruct the reionization history and determine, espe- 
cially, if the reionization process was either instantaneous 
or lengthy and complex as implied by current CMB data. 
Moreover, since observable effects in CMB data depend 
only on the ionized content while 21 cm fluctuations come 
from inhomogeneities in the neutral distribution, the com- 
bination may allow additional information in the study of 
reionization. 

Although there is strong scientific motivation to study 
21 cm fiuctuations during the era of reionization, there 
are several difficulties to overcome. In addition to ex- 
perimental challenges, involving observations at low radio 
frequencies, the background itself is highly contaminated 
by foreground radio emission. Synchrotron and free-free 
emission from the Milky Way (Shaver et al. 1999), low fre- 
quency radio point sources (Di Matteo et al. 2002, 2004) 
and free- free emission from free electrons in the intergalac- 

^http: / /web. haystack. mit.edu/arrays/MWA/site 
^http:/ /www. lofar.org 
■'http:/ /www. skatelescope.org 
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tic medium (Oh 1999; Cooray & Furlanetto 2004) are now 
thought to be the chief sources of confusion. While it has 
been suggested (Oh & Mack 2003) that background stud- 
ies may be impossible due to the large number of contam- 
inating sources (with brightness fluctuations much higher 
than expected for neutral Hydrogen at high redshifts), oth- 
ers have contended that a multifrequency analysis of the 
21 cm data can be used to reduce the foreground con- 
tamination given the smoothness of these contaminants 
in frequency space (Shaver ct al. 1999; Di Mattco ct al. 
2002; Zaldarriaga et al. 2004). Initial calculations suggest 
that, for experiments such as LOFAR and SKA, one can 
easily remove the foregrounds to a sufficiently low level 
that statistical studies related to 21 cm fluctuations can 
be achieved (Zaldarriaga et al. 2004). 

In this paper, we extend the initial discussions related 
to foreground removal by considering a detailed analysis 
of the multifrequency removal technique. When estimat- 
ing how accurately cosmological parameters can be mea- 
sured, we take into account the complication that most 
foreground properties, like the level of smoothness across 
frequency, must be estimated at the same time as the sig- 
nal. Our foreground model includes not only the normal- 
ization, frequency dependence and scale dependence for 
each physical quantity, but also variations in the frequency 
coherence and even the functional form of the foreground 
correlations across frequency. We also consider the extent 
to which 21 cm fluctuations can be cross-correlated be- 
tween frequency bins. This cross-correlation of the signal 
between channels reduces the efficiency with which fore- 
grounds can be cleaned. 

The paper is organized as follows; in §2, we discuss the 
21 cm signal anisotropy and frequency correlations. In 
§3, we give an overview of the most important foreground 
contaminants. In §4 we discuss the experimental setup as- 
sumed in our analysis. In §5 we discuss the method used 
for the error forecast and the multifrequency cleaning tech- 
niques to reduce the foreground contamination. Finally in 
§6 we study how well physical parameters related to the 
21 cm signal can be extracted. Throughout the paper, 
we make use of the WMAP-favored ACDM cosmological 
model (Spergel et al. 2003). 

2. THE 21 CM SIGNAL 

When traveling through a patch of neutral hydrogen, 
the intensity of the CMB radiation will change due to ab- 
sorption and emission. The corresponding change in the 
brightness temperature, T2i(i/), as compared to the CMB 
at an observed frequency v is then 

W^^^^r (1) 

where Ts is the temperature of the source (the spin tem- 
perature of the IGM), z is the redshift corresponding to the 
frequency of observation {1 + z — V2i/i^, with 1^21 = 1420 
MHz) and Tcmb = 2.12,{l + z)K is the CMB temperature 
at redshift z. The optical depth, r, of this patch in the 
hyperfine transition (Field 1959) is given in the limit of 
UbTs » hv2i by 
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where Aiq is the spontaneous emission coefficient for the 
transition (2.85 x 10~^^ s~^) and tihi is the neutral hydro- 
gen density. This can be expressed as rini = a;jjn6(l -|- 5b), 
when nf, is the mean number density of cosmic baryons, 
with a spatially varying ovcrdcnsity Sj, and xh is the frac- 
tion of neutral hydrogen [xh = 1 — -T',; where is the 
fraction of free electrons). We refer the reader to Zaldar- 
riaga et al. (2004) for further details. 

2.1. 21cm, Power Spectrum 

The measured brightness temperature corresponds to a 

convolution of the intrinsic brightness with some response 
function Wv{r) that characterizes the bandwidth of the 
experiment: 



T^(n,z/o)= j 



= / rfrI^,„(r)T2i(n,r), 



(3) 



where n is the direction of observation and r is the ra- 
dial distance corresponding to the observed frequency v. 
We can now determine the angular power spectrum of 
T^{n,VQ) which is related to the 3-d power spectrum of 
T2i(n, r), defined by 

(T2l(k,I/i)T2l(k',Z/2)) = (27r)35^(k + k')P2l(fc,!/l,i/2), 

(4) 

where T2i(n,r) = T2i(r,t/) = / ^T2i(k, z/)e''^"-. Mak- 

ing use of the spherical harmonic moment of the 21 cm 
fluctuations, at a frequency vq, 

atM= JdnYil{h)T%n,uo), (5) 

we write the angular power spectrum as (e.g. Kaiser 1992) 
{atMaZi'^2)} = C!{i.i,U2) (6) 

= ^ y" fc^dfc P21 (fc, z/i , iy2)Ip {k)Ip {k) , 

where 

ir{k)= J drW^ir)jiikr) , (7) 

with the spherical Bessc;! function given by ji{kr). In de- 
riving this form for the angular power spectrum, we have 
made use of the Rayleigh expansion for the plane wave 
given by 



47r^i'iKfcr-)yr(fi)^r(k). 



(8) 



In Figure 1, we show in the top panel the integral, Ij'ik), 
when I = 100 as a fimction of fc, for three different fre- 
quency bins of width IMHz. As the frequency separation 
increases, the oscillations are out of phase suggesting that 
Cf (i^i, z/2) is not significant for 1^2—1^1 ^ lOMHz. However, 
for smaller separations Cf (z^i, 1^2) certainly cannot be con- 
sidered to be zero. This is clear by looking at the bottom 
panel of Figure 1 where we show the signal correlations 
across frequency for separations from 1 to lOMHz. 

To proceed further we need an expression for the power 
spectrum of the 21 cm fluctuations, P2i{k,z), or in this 
case, spatial inhomogeneities in the underlying gas density 
field. In order to simplify the calculation related to fluctu- 
ations in the density field, we assume a scenario where the 
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Fig. 1. — The frequency cross-correlations of the 21cm signal. 
Top panel: the absolute value of the function I'{(t:) when / = 100 
for three frequency bins with a width of IMHz. For large frequency 
separations these functions are out of phase and the integral over k is 
essentially zero. Bottom panel: signal correlations across frequency. 
Each curve corresponds to a frequency separation from IMHz to 
lOMHz (IMHz bandwidth). 



gas density field is heated well above the CMB tempera- 
ture during the reionization process. In addition to reion- 
ization, we also expect the strticture formation process, 
such as shocks in forming virialized dark matter halos, 
where first stars eventually form to subsequently reion- 
ize the universe, to heat the gas such that Tg » Tcmb 
(Gnedin & Shaver 2004). The gravitational shocks will 
also heat the gas in overdense regions like sheets and fil- 
aments such that the assumption we make here with re- 
gards to the gas temperature should be valid throughout 
the IGM at a time around the era of reionization where 
one also expects a substantial background of X-rays from 
supernovae, etc (Venkatesan et al. 2001; Chen & Miralda- 
Escude 2004) to also heat the IGM quickly. At sufficiently 
high z, prior to significant reionization and heating, \vq ex- 
pect Ts < Tcmb and it may be possible to sec HI in 21cm 
absorption (Madau et al. 1997). 

With Ts » Tcmb, the 21 cm signal will be observed 
in emission with respect to the background CMB. Writing 



T2liz) = c{l + (5)(1 - Xe- XeSx), 



(9) 



where 
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and 6x is the perturbation in the ionization fraction {Sx = 
^"g"^" ), the corresponding three-dimensional power spec- 
trum is 

P2l{k, Z) = c2 [(1 - XefPss{k, z) + xlPg^s^ (fc, z) 

- 2PssAk,z)Xe{l-Xe)]. (11) 

The dark matter power spectrum is represented by Pss{k, z) 
(and we are assuming that (5 = (5b), Ps^s^ {k, z) is the power 
spectrum from the perturbations in the ionized fraction 
and Pss^ {k, z) is the cross-correlation power. 

The correlations in the ionized fraction depend on the 
reionization history. To proceed further we used the model 
of patchy reionization in Santos et al. (2003a) to write 



Ps^sAk,z)=b'iz)Ps5ik,z)e- 



(12) 



where b{z) is a mean bias (halo bias weighted by the dif- 
ferent halo properties) and R the mean radius of the HII 
patches. 

We crudely model the time-dependence of R as R = 
( \ 

( 1 Rp where i?p is the (comoving) size of the fun- 
damental patch which wc take to be i?p ^ 100 Kpc. The 
cutoff scale R is equal to i?p when reionization starts and 
increases with time as HII regions overlap and form larger 
HII regions. The dependence of R/Rp on Xe follows from 
assuming a Poisson distribution of patches and that their 
volumes add when they overlap. The essential consequence 
of this dependence of R on x^ and our choice of Rp is that 
for I < 10, 000 (scales greater than 1 Mpc) the transition 
is very sudden; i.e., R increases from 1 Mpc to oo very 
rapidly. 

In reality there is a distribution of patch sizes and there 
are spatial correlations between ionizing sources which cause 
some of the patches to overlap earlier than they would 
given a Poisson distribution. Furlanetto et al. (2004b) 
have recently developed a sophisticated model of the size 
distribution of patches and their correlations. They find 
a much more gradual transition to large patch sizes. In 
their model, large patch sizes suppress fluctuation power 
even as low as / = 1000 while the ionization fraction is still 
as low as Xe = 0.5 (although this depends on the assumed 
efficiency with which the collapsed baryons ionize the sur- 
rounding hydrogen). We caution the reader that modeling 
the signal on scales smaller than the largest patches is dif- 
ficult and the shape of the signal on these scales remains 
highly uncertain^. 

The bias in equation 12 is much larger than 1 so that 
wc can safely neglect the cross-correlation contribution in 
our analysis. We show in Figure 2 the 21cm power spectra 
for three reionization models with the same optical depth 
(r = 0.17). Typically, for ^; > 7 the reionization models 
give R ^ 0.2Mpc so that the size of the patches are only 
important for I > 10^ and below that, the power spectrum 

^Note also that for a IMHz bandwidth the signal will be smoothed 
out on scales smaller than ~ 18Mpc at z = 10. 
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Fig. 2. — The power spectrum for the 21cm signal at z=9.2 and 
assuming a bandwidth of Ai^ = IMHz. The solid, dashed and dot- 
dashed lines use the corresponding reionization models shown in 
Figures 3 and 4. 
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Fig. 3. — Xe for three different fiducial models with r = 0.17. 
Solid line is for the fiducial model used throughout the paper. 
Dashed line assumes that H2 feedback effectively shuts off star for- 
mation in all Type 11 halos. Dot-dashed line assumes that a transi- 
tion from metal free to normal star formation takes place abruptly 
at a certain redshift. 

just follows the shape of the matter density perturbations. 
Note that the bend in the curves at I ^ r/6r ^ 500, where 
6r is the width of the window function, Wi,{r), is due 
to the Limber regime where perturbations are smoothed 
out. Figures 3 and 4 show the ionization fraction and bias 
for the corresponding reionization models from Haiman & 
Holder (2003). 



Fig. 4. — Effective bias for the same fiducial models as in Figure 3. 



3. FOREGROUNDS 

We consider four different foregrounds: Galactic syn- 
chrotron, Galactic free-free, extra-galactic diffuse free-free 
and extra-galactic point sources. All of these contaminants 
produce much more fluctuation power at each frequency 
than the 21cm signal. Measurement of the 21cm signal 
would be impossible if not for the high coherence of the 
contaminants across frequency, compared to the very short 
frequency space correlation length of the signal. 

With i and j labeling the different contaminants let us 
define 

Ci^{u,,i^2)^{aiMaiM. (13) 

We expect correlations between the different contaminants 
to be negligible. We model the foregrounds as power laws 
in both I and v so that: 

Cf{u,i^) = Ai{1000/lf'{uf/iyf^'. (14) 

where = 130MHz. We also expect the foregrounds to 
be highly coherent, i.e. I"{i>i, 1^2) — 1 where 

ir{vuV2) = Cr{v^,V2)/^Cf{vuVl)Cf{v2,V2). (15) 

Althotigh we consider other parameterizations of our 

starting point is 

ir{u^,V2) = exp [- log2(i.i/z.2)/2/a] , (16) 

which, for the frequency range and values of we will be 
considering, reduces to 

/f(^/i,i^2)^i-iog'(j^i/;^2)/2/ei;. (17) 

This form describes the departure from complete correla- 
tion to lowest non-trivial order. It follows if one assumes 
the underlying sources have power-law spectra with vary- 
ing spectral indices a. If sources with different spectral in- 
dices have C;'s with different shapes (for example, because 
they have different redshift distributions), then /; can be 



5 



/-dependent. Note that, due to this decorrelation, the fre- 
quency dependence of the foreground brightness temper- 
ature will actually have departures from a power law and 
change with position on the sky. 

If Ci is dominated by Poisson fluctuations, then ^ = 
1/Aa where Aa = ((a— a)^)-'^/^. This was the case consid- 
ered by Tegmark (1998). If C/ is dominated by clustering 
and d In C//c? In a — then ^ — Gaj ^o? where, following 
Zaldarriaga et al. (2004), we assume that correlations be- 
tween sources with a\ and a-2, fall off as exp \—[a\ — a-i)^ j (2(t^ 

In Table 1 we show the parameters of our fiducial models 
for the four different contaminants that we now discuss. 
Figure 5 shows the expected foreground power spectrum 
together with the signal at 140MHz. 

3.1. Extragalactic point sources 

The extragalactic sources of contamination likely to be 
important at the relevant low flux levels are radio galaxies, 
active galactic nuclei and normal galaxies. There are sev- 
eral surveys of extragalactic sources in or near the wave- 
length range of interest including the 6C survey (Hales 
et al. 1988) at 151 MHz, the SCR survey at 178 MHz 
(Laing et al. 1983) and a survey by Cohen et al. (2004) at 
74 MHz. 

All but the Cohen et al. (2004) survey were used by Di 
Matteo et al. (2002) to estimate the angular power spectra 
from unresolved sources at 150 MHz. They model both the 
Poisson contribution and the correlated contribution, as- 
suming the sources are clustered like Lyman break galax- 
ies (Giavalisco et al. 1998) at z ~ 3. They find that if 
sources brighter than S'cut = 0.1 Jy can be detected and 
removed, then the clustering signal, and not the Poisson 
signal, will be the dominant contribution to the variance of 
maps smoothed on angular scales larger than about 1'. As 
the flux cut is reduced further, the Poisson contribution 
drops faster than the clustering contribution, due to the 
Poisson contribution's heavier weighting of the brightest 
sources. 

We use the Di Matteo et al. (2002) model with S'cut = 
0.1 mjy for the angular power spectrum of the background 
of unresolved sources ^. This translates into A = 57mK^, 
and /3 = 1.1 as shown in Table 1. 

For the frequency dependence we turn to Cohen et al. 
(2004) who compare fluxes of 947 of their 949 sources that 
appear also in the NVSS survey at 1400 MHz (Condon 
et al. 1998). Assuming the spectrum of each source follows 
a power-law in intensity, Ii, oc v'^, they find an average 
power-law index of 774"° = —0.79 with a scatter of about 

^Note that the model for the differential source counts is dN/dS = 
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extragalactic point sources 


57.0 


1.1 


2.07 


1.0 


extragalactic free-free 


0.014 


1.0 


2.10 


35 


Galactic synchrotron 


700 


2.4 


2.80 


4.0 


Galactic free-free 


0.088 


3.0 


2.15 


35 



Table 1 

Fiducial parameters at I = 1000 and 130MHz. 
Parameters are defined in equations 14, 15 and 16. 



(57 = 0.2. They also see a trend for flattening of spectral 
index as the 74 MHz flux density drops from 1 Jy down to 
0.2 Jy going roughly as 

^ifo = -0.83 - 0.131ogio (174/ Jy) . (18) 

Extrapolating further down to 0.1 mJy we find 774"° = 
-0.31. 

Cohen et al. (2004) also find 545 of the 947 sources in the 
WENSS catalog at 325 MHz (Rengehnk et al. 1997). They 
find a significant flattening of the spectra of most sources 
towards lower frequencies. The mean spectral indices are 
related by 7|f = 0.24 + 7^|™ = -0.07. 

Since fluctuations in antenna temperature follow as 5T{v) oc 
v~'^5Iy, then a = 2 — 7. Our best guess based on the ex- 
trapolations to lower flux densities and lower frequencies 
is that 7 = —0.07 and therefore we use a = 2.07 as our 
fiducial value. 

Since clustering, according to Di Matteo et al. (2002), 
is the dominant source of fluctuation power we have ^ = 
(Ta/Aa^. For our fiducial model we take Aa ^ 5^ ~ 0.2, 
consistent with Cohen et al. (2004). The value for aa is less 
clear. As Zaldarriaga et al. (2004) argue, if the sources are 
all tracing the same dark matter distribution, even with 
different biases, then they would be perfectly correlated 
with each other; i.e., CTq = 00. Certainly on large scales we 
expect this to be a good approximation. In the context of 
the halo model (Cooray & Sheth 2002), the relevant length 
scale is the size of a typical halo hosting these sources. How 
this translates to an angular scale depends on distance to 
the source. 

We have not attempted to calculate ctq here, but merely 
point out that even taking the conservative assumption 
that ctq, = Aa we get ^ = 5. Because, as we will see, we 
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Fig. 5. — The 21cm power spectrum (blue thick line) at z=9.2 
{v = 140 MHz) and its foregrounds. Point sources - red thin solid 
line, Galactic Synchrotron - red thin dotted line, Extragalactic free- 
free - red thin dashed line and Galactic free-free - red thin long- 
dashed line. CMB - green dot-dashed line. 
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can afford to be even more conservative, we set our fiducial 
value of ^ even smaller to ^ = 1. 

We note that although for low S'cut clustering may be 
the dominant source of fluctuation power, Poisson fluctua- 
tions may remain as the dominant source of dccorrclation, 
especially if cTa >> Aa. In this case though, ^ would still 
be greater than 1/Aa and therefore even more coherent 
than in our fiducial model. 

We also note that although the extrapolations from cur- 
rent observations of sources in the relevant wavelengths are 
over many decades in intensity, we are, of course, not stuck 
with these extrapolations. With the type of experiments 
we consider here, we will be able to study large numbers 
of resolved sources with flux densities just above S'cut- 

3.2. Extragalactic free-free emission 

The low-frequency radio background produced by ion- 
ized electrons, via the free-free emission, is known to be 
a significant foreground for 21 cm studies (Oh 1999; Oh 
& Mack 2003; Cooray & Furlanetto 2004). The free-free 
emission coefficient can be written as 

= 5.4 X 10-39n2T-i/2g//(:.,T,)e-''^/'=^' 

ergs cm~^ s"^ Hz"^ sr~^ , (19) 

where the Gaunt factor can be approximated in the radio 

regime as .g//(t/,Te) « 11.96Tf i^i^-o i (Lang 1999). The 
cumulative specific intensity is simply = J dxi-v{z)x^{z) 
C{z)/{1 + zY, where x is the conformal time and v{z) = 
^'obs(l + z) with z/obs the observed frcqiiency. The bright- 
ness temperature is = c?I,j2v'^T^Jf,B and the electron 
fraction Xc{z) captures the mean ionization history of the 
IGM. The clumping factor of the electron density field is 
defined as C{z) = (ng)/(ne)^. When estimating the free- 
free background and its spatial fluctuations, this is the 
most uncertain parameter. Initial estimates of the free-free 
background suggest that gas clumping boosts the specific 
intensity by a factor from order unity at 2; ~ 20 to over 
100 at z < 3 (Oh 1999). 

A variety of models from the literature considered in 
Cooray & Furlanetto (2004) suggest that the free- free back- 
ground is uncertain by at least 2 orders of magnitude. 
Among these models, a favorite is to estimate the clumping 
from the large-scale matter distribution by, for example, 
assuming that ionized electrons reside in virializcd dark 
matter halos (Haiman et al. 2001; Benson et al. 2001; 
Gnedin & Ostriker 1997) though such models are only 
expected to yield a minimal estimate for C{z), because 
the models do not include recombinations within galax- 
ies (where most ionizing photons are probably absorbed) 
or the biased distribution of ionizing sources. Numerical 
simulations by Gnedin (2000) show that these two fac- 
tors strongly affect the effective clumping factor early in 
the reionization process, when ionizing photons are con- 
fined to the dense regions around galaxies. The second 
approach, taken by Oh (1999), is to construct a model for 
the production rate of ionizing photons nion and assume 
that, throughout the universe, ionization equilibrium is a 
good approximation. The quantity x1C{z) is then fixed 
through the relation 

hion = xl{z)nl{z)C{z)arec, (20) 
where arec is the recombination coefficient. This approach 
predicts stronger clumping than models based on the IGM 



gas as it includes recombinations inside galaxies. These 
dominate the clumping so long as the escape fraction of 
ionizing photons is small. However, this method depends 
on the assumed (and uncertain) ionizing photon produc- 
tion rate. It also docs not include coUisional ionization in 
massive halos at low redshift, though this is expected to 
be small. 

The results related to these two approaches and a com- 
parison to results from numerical estimates of gas clump- 
ing are summarized in Cooray & Furlanetto (2004). In 
Fig. 5, we make use of the high-end estimate of free-free 
anisotropy spectrum based on a model related to ioniza- 
tion equilibrium and using the star-formation history de- 
rived by Hernquist & Springel (2003). In addition to the 
angular power spectrum, the foreground model requires in- 
formation related to the coherence of fiuctuations between 
frequencies. Given that the extragalactic free- free emis- 
sion is dominated by point sources at z ^ 3 which follow a 
biased distribution with respect to the linear density field, 
the angular power spectrum is dominated by the cluster- 
ing nature of this point source distribution. The param- 
eter of interest, ^, then takes the value of Uaj I^o? where 
a is the frequency dependence determined by a combina- 
tion of the Gaunt factor and Q-hv/kT^_ Even accounting 
for variations in the electron temperature by an order of 
magnitude, from lO'' to 10^ K, the spectral index, as a 
function of frequency, varies slowly, since hv <C kT^. We 
take this small variations into account by sotting ^ = 35 
(with Aa ~ 0.03 assuming cxa « Aa). This value is sig- 
nificantly high when compared to, say, radio point sources 
with a value of ^ = 1 in our fiducial model. The difference 
is due to the fact that observed spectral indices of point 
sources are rapidly varying from one source to another 
while free-free frequency spectra, from source to source, is 
expected to be very uniform due to the nature of free-free 
emission. 

While the anisotropy spectrum and frequency coherence 
are just estimates based on model calculations, unfortu- 
nately, there are no observational measurements of the 
free-free background either in terms of mean brightness 
temperature or spatial fluctuations. The planned Abso- 
lute Radiometer for Cosmology, Astrophysics, and Diffuse 
Emission (ARCADE) experiment may be able to study 
some aspects of free- free background at frequencies around 
a few GHz. Also, low- frequency experiments such as SKA 
will have a high frcriuency imaging capabilities such that 
one can use data at frequencies just above 1420 MHz to 
determine the extent to which free-free emission can be a 
contaminant at low radio frequencies. While the free- free 
background can be a problem, in Cooray & Furlanetto 
(2004) it was shown that the integrated signal is domi- 
nated by low- 2; sources (at redshifts less than 3). Also, a 
large fraction of the free-free fluctuations comes from rel- 
atively bright point sources that can be cleaned from the 
maps. Given that most of the emission is at low redshifts, 
when compared to 21 cm fluctuations from the reioniza- 
tion era, the free-free spectrum is smooth in frequency. 
Thus, multi-frequency differencing is expected to remove 
the smooth background radiation but not 21 cm radiation 
as the two fields are strongly uncorrelated given the dis- 
joint in redshift distributions; we do expect some free-free 
emission from the era of reionization, but this was shown 
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to be substantially below the noise level of upcoming ex- 
periments in Cooray & Furlanetto (2004). 

If the fraction of emission at ^; ~ 10 were to be higher 
for some reason, then the confusion with free-free fluctua- 
tions may become substantial; for free-free fluctuations to 
dominate 21 cm anisotropics, even with optimal multifre- 
quency cleaning, we require the fractional contribution to 
total emission from 2 ~ 10 to be at the same level as 2 ~ 1, 
which we do not consider to be physically possible. Thus, 
it is unlikely that after multifrcquency cleaning, free-free 
fluctuations will be the dominant residual contamination 
of 21 cm maps. 

3.3. Galactic Backgrounds 

While radio point sources and ionized halos that emit 

free- free emission are expected to produce small-scale fore- 
ground anisotropy structure at low radio frequencies, the 
large angular scale fluctuations are expected to be dom- 
inated by the synchrotron emission within the Galaxy. 
The synchrotron background is moderately well under- 
stood, when compared to say low frequency faint radio 
sources or the free-free background, given that it also can 
affect temperature anisotropy measurements of the CMB 
at frequencies around 30 GHz. In addition to temperature 
anisotropics, the highly sensitive CMB polarization ob- 
servations, with substantially smaller signals, also require 
a thorough understanding of the synchrotron background 
from the Galaxy. 

In Fig. 5, we make use of a "middle of the road" es- 
timate for the Galactic synchrotron fluctuations, follow- 
ing Tegmark et al. (2000). The synchrotron spectrum 
has an angular power spectrum that scales with multipole 
as Ci ~ Z~^ '* (Tegmark et al. 2000), though this slope 
is rather uncertain with observational measurements, e.g. 
from the 408 MHz Haslam map, varying from -2.5 to - 
3.0 down to angular scales of a degree. In addition to 
the spatial spectrum, we take a frequency spectrum that 
scales as h'^" and take a = 2.8 following the data from 
Platania et al. (1998). The frequency spectrum, however, 
varies across the sky with Aa ~ 0.15 (Tegmark et al. 
2000), though information related to aa is not yet avail- 
able from the data. We take a conservative estimate of 0.1 
and set ^ ~ Uaj Ao?, as large scale clustering dominates 
the anisotropy spectrum, to be 4 in the case of Galactic 
synchrotron fluctuations. 

Note that the synchrotron amplitude is set based on 
foreground analysis of CMB experimental data at 19 GHz 
(cross-correlated with low frequency maps) where the power 
variance of synchrotron fluctuations was estimated to be 
52± 17 /iK at an angular scale of 3 degrees by de Oliveira- 
Costa et al. (1998). Though fluctuation power decreases 
with increasing angular resolution, the amplitude for these 
fluctuations is such that the synchrotron emission from the 
Galaxy dominates temperature fluctuations in all angular 
scales of interest. 

Fortunately, the synchrotron emission can also be con- 
tained in upc;oming 21 cm observations given the fact that 
this power is expected to be smooth over a wide band in 
frequency space while 21 cm fluctuations related to neu- 
tral Hydrogen at redshifts around reionization will vary 
significantly. 

In addition to synchrotron emission, the Galactic free- 



free background is also expected to confuse low-frequency 
radio measurements of the neutral Hydrogen distribution 
around the era of reionization. Following Tegmark et al. 
(2000), we include this contribution as a separate compo- 
nent with a power spectrum that scales as C; ^ and 
a frequency spectrum with v~'^-^^ . As in the case of extra- 
galactic free-free background, we assume a large frequency 
coherence and set ^ = 35 when Aa = 0.03. Note that the 
Galactic free-free background has a similar spatial power 
spe(;trum as that of the synchrotron background, but with 
a lower amplitude. While such a background can be ig- 
nored, given the lower amplitude compared to Galactic 
synchrotron, we include this component due to the differ- 
ence in frequency dependence and the fact that it has the 
same frequency scaling as that of the extragalactic free- free 
component. By doing so, we can investigate how differ- 
ences in frequency scaling and the power spectrum scaling 
can complicate a potential frequency cleaning attempt. 

4. EXPERIMENTAL SETUP 

For an interferometer type experiment, the measured 
flux in a given visibility ^(u, z/), can be expressed as 

y(u, y) = ^-§^\ dnT{n, ,y)A{n, i/)e2-"-", (21) 

where u corresponds to a point in the u — v plane, T(n) is 
the sky temperature, A{n) is the primary beam and 
converts temperature to flux (see e.g. White et al. 1999). 
The visibilities will then consist of a convolution of the 
Fourier transform of T{n) 

y(u, y) = ^J d?u'f{n', v)A{vi - u', v) (22) 

with T(u) and A(u) the Fourier transforms of T(n) and 
A(n) respectively. The typical width of ^(u, u) will be set 
by the size of the dishes d, so that Am = d/X and visibil- 
ities closer than this in the w, — v plane will be strongly 
correlated. The Fourier component T(u) is basically the 
flat-sky version of the spherical harmonic transform of the 
sky temperature, aim (for the correspondence see e.g. San- 
tos et al. 2003b). Its power spectrum is given by 

(f (u,i/)f*(u',z/')) = ^^(u- u')C;=2^|u|(i^,i'')- (23) 
Finally, we can relate the power in the visibilities with 
the usual angular power spectrum (with the relation / = 
27r|u|): 

X J d=^u'i(u-u',I/i)i*(u-u',Z/2)Ci=2,r|u'|('^l.'^2) 

~<-'(=2^iu|(i'i,j^2) j ) 

X j d^u'A{vi-u',ui)A*{n-M',V2), (24) 

where the last line is a reasonable approximation since the 
power spectrum Ci is quite smooth. From this expres- 
sion we can see that, because the beam size depends on 
the frequency, the visibilities will be less correlated across 
frequency than the underlying signal on the sky (Di Mat- 
teo et al. 2002). Basically the width of A{vl) increases 
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with frequency so that visibihties at higher frequencies will 
have contributions from a larger range of T(u) (equation 
22). This is particularly important in the foreground case, 
since their correlation structure is crucial for the cleaning 



process. Using / cPu A{u,iy) 



J£u' |A(u-u',i/)|2 ^ 
tions can be written as 



{V"iu,i^i)V' 



l/{Au) 

(U,i^2)) 



1 - (Au)2A(0) so that 
, the foreground correla- 



1^1 



'i=2Tr|u| 



(1^1,1^2)-, (25) 



1^2 



V(|y»(u,I.l)|2)(|F"(u,i.2)P) 

where I"{vi,V2) is given by equation (15) and 1^2 > J^i- 
Even if (1^1,1^2) = 1 the correlation can still be much 
smaller than 1 , making foreground removal extremely dif- 
ficult. It will therefore be essential for this type of ex- 
periment to generate an effective beam size that is con- 
stant across frequency. These considerations will place 
stringent constraints on experimental design and observ- 
ing strategy. Only with sufficiently dense sampling and 
with the primary beam known sufficiently accurately can 
higher frequency visibilities be deconvolved to match the 
lower frequency visibilities. From now on we will assume 
this is possible and drop the explicit frequency dependence 
of the beam. 

The corresponding noise in each visibility is given by 
(Rohlfs & Wilson 2004) 

2 



(|7V(u)p) = 



sys 



A. 



dish 



1 



(26) 



where T^ys is the system temperature, Az/ is the band- 
width, and Adish is the area of each individual antenna in 
the array. From equation (24), the noise angular power 
spectrum will then be 
2 



1 



u' 



(27) 



Adish ) Aiytv J d'^u' \A{u- 

Moreover, assuming the baselines are arranged so that the 
Fourier coverage is uniform, the time spent for each vis- 
ibility will be roughly constant. At a given instant, the 
experiment will observe an area in Fourier space given by 
Nb{Au)'^, where Nb is the number of baselines. The ob- 
servation time per unit area will then be 

= ^" (28) 

iVb(Aw)2 TTM^^,,^' 

where Umax corresponds to the largest visibility observed 
and to is the total time spent observing a given patch of 
the sky. This way, using Nb ~ 



^sys\ 



where fco ver 
the dishes: 

fct 



A , j-2 J2 ' i'^^) 
^^^^J cover^max 

is the fraction of the total area covered by 



NdishAdish _ Ndishi^i)'^ 



A 



total 



max 



(30) 



For our analysis, we used fcover — 0.034, Imax — 5000, 
Tsys = 200K and Av = IMHz. We choose the total ob- 
servation time to be ttotai = 4 months. To specify we 
took into account that it might be more efficient to observe 
different patches of the sky in order to reduce sample vari- 
ance. The optimal observation time for one patch, given 
the above values, is to '^-^ 1 month. Note that by optimiz- 
ing the signal to noise (signal/noise ~ 1), the expected 



errors will scale as ~ 1/ fcover/ Vttotai, so that it should 
be straightforward to extrapolate our results to other ex- 
periments (see Morales 2005 for a discussion on the design 
of Epoch of Reionization observatories). The correlation 
length across multipoles is set by the dish size, which we 
assumed to be d ~ 33 m corresponding in multipole space 
to AZ = 27rAu ^ 100 at 2: = 9. Figure 6 shows the noise 
level we use, together with the expected signal at 140MHz. 



5. ERROR FORECAST 

Our goal is now to calculate the expected error in the 
21cm power spectrum, taking into account the experimen- 
tal noise and foreground contamination. Given a set of 
measurements {y}, the information of how well we can 
measure a set of parameters {p}, is contained in the Fisher 
information matrix (Kendall & Stuart 1979) 

-'-my 

where C{y,p) is the likelihood of parameters {p} given 
observations {y} and the matrix is to be evaluated at 
the true values of the parameters {p}. For now, we con- 
sider these parameters to be the 21cm auto power spectra 
binned in multipole and frequency space plus the fore- 
ground model parameters. The Cramer-Rao inequality 
states that {F~^)ij is the smallest variance that any unbi- 
ased estimator of the parameters can have (see Tegmark 
et al. 1997 for a review). 

For an interferometer, we can take the set of measure- 
ments {y} to be the measured visibilities, V{u, 1/). Assum- 
ing all fluctuations are Gaussian with zero mean ((y) = 0) 
the Likelihood is entirely specified by the covariance ma- 
trix C = (yy*) = {V{u,iy)V*{u',iy')). Note that the 
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Fig. 6. — The noise power spectrum assumed in our analysis (red 
dashed hne). Blue solid line - the 21cm signal at z=9.2 (140 MHz) 
with Au = 1 MHz. 
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brightness temperature is not exactly Gaussian (see Furlan- 
etto et al. 2004a). In order to obtain large values of the 
temperature we need regions with large densities, but these 
are typically ionized giving T21 — 0, so that the Gaus- 
sian distribution will be cut for large values. However, as 
long as we are probing scales much larger than the effec- 
tive patch size, the signal will be an average over many 
bubbles and we expect the Gaussian assumption to be a 
good approximation for the Fisher matrix analysis. We 
also make the approximation that the foreground fluctu- 
ations are Gaussian distributed. While this assumption 
holds true for high-redshift contributions, such as free-free 
point sources, foregrounds such as synchrotron emission 
within the Galaxy are highly non-Gaussian over large ar- 
eas on the sky. Our Gaussian assumption is then expected 
to underestimate the effect of foregrounds at large angular 
scales corresponding to tens of degrees. However, here we 
will consider measurements of 21 cm fluctuations at a de- 
gree scale and below so that the Gaussian approximation 
is expected to be a reasonable one. 

The covariance matrix can be quite complex due to the 
cross-correlations between different u. To proceed further 
we assumed a bin size in I = 2it\u\ space such that this 
cross-correlation is negligible. The matrix C will then be- 
come block diagonal, making the Fisher matrix analysis 
more straightforward. As discussed in the previous sec- 
tion, it will be necessary to create an effective beam such 
that the correlation length is constant across frequency. 
We choose Al = 100 which, for the experimental setup de- 
scribed above (§4) makes the values almost uncorrelated 
among different I. The number of independent measure- 
ments at a given I = 27r|u| is then Ni = ^Aito"' • 

We will proceed the analysis making reference to the an- 
gular power spectrum instead of the power in the visibili- 
ties, knowing that the two can be related through equation 
(24). The covariance at a given I between two frequency 
bins is: 

Cl^ = Ciiui, u^) + Ctiy,, v^) + Cf{v,, y,)6,, . (32) 
The corresponding Fisher matrix is then. 



(c')-^|^(c')-^|^ 

dPa dpb 



(33) 



where {pa, Pb} are the parameters we wish to estimate. 

5.1. Power Spectrum 

In this section we investigate how well we will be able 
to measure the 21cm power spectrum. To better under- 
stand the impact our foreground models will have on fore- 
ground removal, we start by just considering one type of 
foreground (point sources) and assume that all its statisti- 
cal properties are known. Our parameters are then: pa = 
Ci^iyai ^a) and we assume that the correlation of this sig- 
nal across frequency, c. g. Cf {Va , Vb)l\JCl (Va , Va)C^ {vb, l^b) , 
is known from theory. 

We can calculate an approximate expression for the ex- 
pected error in Cf^{iyi, Ui) which will make it easier to un- 
derstand the results from the Fisher matrix analysis. If the 
foregrounds were perfectly correlated across frequencies 
{ci{^(v2) = I'^imi'^i)) ^^'^ signal uncorrelated, then 
clearly by subtracting the two maps (a;m(f2) — jaim{t^i)) 
the foregrounds would be completely removed. A decrease 



in the foreground correlation translates into fluctuations 
in 7 so that the foreground subtraction is not as effective. 
Using only two frequency bins, the corresponding error of 
the 21cm power spectrum is then (see Zaldarriaga et al. 
2004) 

2 r , -,2 



21\2 



[Acn 



Cf 



I) 



(34) 



2^ + ir'+^'" + ^Ai 

where / is the correlation coefficient and we have assumed 
the signal power spectrum to be the same in both frequen- 
cies. Cleaning in this case is very effective, with foreground 
contamination down by Cf {1 — I). 

We want to reconstruct the signal power spectrum as 
a function of frequency and so now drop this assump- 
tion of frequency independence. Without this assump- 
tion, two channels arc now only sufficient for constraining 
a linear combination of the signal power spectra, not the 
signal power spectra themselves. However, with N chan- 
nels, N{N — l)/2 different foreground- cleaned maps can be 
made (one for each pair of maps) which is more than the 
N signal power spectra to be determined if > 2, mak- 
ing the reconstruction possible. The more channels, the 
more pairs of maps and the more estimates of the power 
spectra we have. If all of these map pairs had similar, and 
independent, information about the signal power spectra, 
then we would expect the errors on each power spectrum to 
improve as y/{N — l)/2. Actually, the errors are corre- 
lated, and the difference maps made from larger frequency 
separations have greater residual foreground contamina- 
tion, so the improvement is slower than (N — l)/2 if 
one increases N by increasing the frequency coverage. 

The growing decoherence of the foregrounds over large 
frequencies limits how well an expanded frequency cover- 
age can improve the reconstruction. Likewise, the coher- 
ence of the signal over very small frequency separations 
limits how well a refined binning can improve the recon- 
struction. With a fine enough binning, moving to finer 
binning no longer increases the effective number of in- 
dependent frequency bins. Figure 7 shows the expected 
foreground and signal correlations across frequency. Al- 
though, as expected, the foreground is much more strongly 
correlated than the signal, there is a non-negligible sig- 
nal correlation. The moderate correlation strength for fre- 
quency separations equal to the bandwidth (about 30%) 
suggests that there will be some improvement in recon- 
struction from making bands finer than 1 MHz, but that 
1 MHz is near the point of highly diminished returns, at 
least for I < 1000. Binning an order of magnitude more 
finely than our nominal 1 MHz bins would also result in 
the added challenge of having to take into account redshift 
distortions due to the peculiar velocities of the neutral gas. 
With 1 MHz binning the redshift distortions are negligible. 

In Figure 8 we show the expected error in the 21 cm 
power spectrum taking into account the experimental setup 
described above and only assuming point source contam- 
ination. We used the point source parameters from Ta- 
ble 1. The red thin lines use information from 40 fre- 
quency bins (from 120MHz to 160MHz) while the green 
thick curves only use ten bins (from 135MHz to 145MHz). 
The larger the number of bins, the smaller is the recov- 
ered error, as expected from the above discussion. Since 
foreground subtraction is diflacult at frequency separations 
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less than IMHz, the best way to improve on this result 
will be to expand the frequency range, although this can 
be problematic due to the man-made electromagnetic in- 
terference. Note however that the recovered error using 
our fiducial point source model (thin solid line) is already 
quite close to the ideal situation with no foreground con- 
tamination (dotted line). Assuming perfectly correlated 
foregrounds (short-dashed line) , the errors come very close 
to the no foreground case; they would be identical were it 
not for the signal correlation. Comparing the solid and 
long-dashed curves, we also see that the Gaussian form of 
the foreground correlations (e.g. equation 16) allows the 
foreground removal to be more efhcient than in the con- 
stant correlation case. This is an important result that we 
will explore in the next section. 

Note also that we are presenting this result for a bin 
at the middle of the frequency range. Figure 9 shows the 
power spectrum and its errors at I — 1000, versus the fre- 
quency, again assuming only point source contamination. 
As we move closer to the end of the frequency interval, 
the error increases, since the number of neighboring chan- 
nels decreases (also, for lower frequencies, the foreground 
contribution is typically larger). However, using the full 
frequency range still provides great statistical power when 
distinguishing between different reionization models as we 
will see later in §6. 

5.2. Dependence on assumptions about frequency 
coherence and its functional form 

In general, the more coherent a foreground is, the easier 
it will be to remove it. This can be seen in Figure 10 where 
we plot the recovered error for the power spectrum at I = 
1000 as a function of the correlation length, ^ (solid curve). 
As above, we are using a Gaussian form for the correlation 
function. As expected, the situation gets worse as we move 
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Fig. 8. — Expected error on the power spectrum at 140MHz 
(IMHz bin size), assuming only point sources contamination (no 
free parameters). Blued dot-dashed line - the 21cm signal. All other 
curves show the recovered errors for different assumptions about the 
point source model, using 40 frequency bins (thin red lines) and 10 
bins (thick green lines). Dotted line - error with no foregrounds. 
Solid curve - assuming a Gaussian correlation across frequency. 
Short-dashed curve - point sources fully correlated. Long-dashed 
curve - assuming a constant correlation. This constant correlation 
case uses / ~ 1 — 3 X 10^^ (the value of the Gaussian at IMHz 
separations). 




Fig. 7. — Expected frequency cross-correlations for foregrounds 
(top) and signal (bottom) at I = 1000 (Ai^ = 1 MHz). 



u 




120 



130 



140 



150 



160 



u (MHz) 



Fig. 9. — The 21cm signal at I = 1000 as a function of frequency 
(blue solid curve) and the corresponding errors. We are assuming 
only point source contamination (parameters from Table 1). Red 
thin dashed line shows the recovered errors using 40 frequency bins. 
Green thick dashed line shows the same for 10 bins. 
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from the ideal case of perfect foreground coherence (large 
^) to completely incoherent foregrounds (small ^). In the 
latter case, the foregrounds simply act as a huge noise 
component. Foreground removal is shown to be effective 
as long as ^ > 0.1 which, if the point sources have a Poisson 
distribution, would translate into Aa < 10. If the point 
source fluctuation power is mainly from correlations and 
we conservatively set CTq, = Aa, then ^ > 0.1 once again 
translates into Aa < 10. We see that our extrapolations 
based on data with brighter point sources would have to 
be very far off in order for the uncertainty in the 21cm 
fluctuation power spectrum to be significantly increased. 

The shape of the correlation function (equation 15) also 
has a very important impact on how well foregrounds can 
be removed. In fact we have already seen in Figure 8 that 
the (unphysical) assumption of constant cross-frequency 
correlations results in larger errors, even though the fore- 
ground correlations between different frequencies are al- 
ways larger than in the Gaussian case (a similar effect, al- 
beit in a different scenario, was already noticed in Tegmark 
et al. 2000). Since the coherence function, I{yi,V2), has 
not yet had its shape accurately measured, we decided to 
repeat the Fisher analysis for a variety of such functions. 
We can rewrite the coherence function (equation 16) as 
7(1/1,1/2) = / \\-0g{v\/v2)/^], where, for our fiducial model, 
f{x) ~ Although this function should behave like 

a parabola near the origin ( equation 17), with /(O) — 1, 
/'(O) = and /"(O) = -1 (Tegmark 1998), it does not 
have to necessarily be a Gaussian, so that we also consid- 
ered functions of the form 




Fig. 10. — Dependence of the expected error on the foreground 
correlation length (at z=9.2 and I = 1000) using a Gaussian cor- 
relation function (blue solid line). Red dotted line shows the same 
for the thick solid curve in Figure 11 (e.g. the worst case scenario). 
Blue dashed line shows the signal. 



We recover the Gaussian case for n = 00 (actually, for the 
typical range of values in x we are considering, n = 10 is 
already a good approximation to a Gaussian). 

Figure 11 shows the foreground correlations for several 
values of n (decreasing n increases the values at the tails 
of the distribution). In Figure 12 we show the recovered 
error in the power spectrum at / = 1000 using the above 
functions. In this case we are only using 10 bins of 1 MHz 
each, centered around 140MHz (solid line). Notice that, 
although all the functions have larger cross-frequency cor- 
relations than the Gaussian case, foreground removal is 
always worse than with the Gaussian form (large n). The 
behavior of f{x) away from the origin is quite important. 
As n decreases, the second derivative of f{x) for x 7^ 
also decreases, so that the function approaches the linear 
(and almost constant) case. This translates into an in- 
crease in the recovered error which peaks at n 3 x 10~^. 
At this point, the correlation for separations of IMHz 
(the bin size we have been using) is still approximately 
the same as in the Gaussian case. As n decreases fur- 
ther, the frequency correlations increase (even at IMHz 
separations) and the function converges to 1 everywhere 
(foregrounds fully correlated), so that the recovered error 
decreases again. Clearly, for foreground removal to be ef- 
fective in the worst case scenario function, we need the 
correlation length to be larger. Figure 10 also shows the 
recovered error as a function of correlation length for this 
worst case value of n (red dotted lines). In order to mea- 
sure the signal we need > 1 which implies Aa < 1. 

5.3. Simultaneous estimation of foregrounds and the 
21cm power spectrum 
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Fig. 11. — Several functions for the foreground correlations of the 
type: (1 + x^/2/n)~"' with x = log(y\lvi) j i^. Gaussian: n^ oo 

(e~^ /2). Red thick solid line gives the largest error and is very 
similar to the constant case with a correlation at IMHz separations 
identical to the Gaussian case. 
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Up to now we have assumed that aU the statistical prop- 
erties of the foregrounds are perfectly known. In fact, we 
hope to learn a lot about the spatial and frequency depen- 
dence of foregrounds from experiments set to measure the 
21cm signal. Information about the spectral index of point 
sources at the frequencies of interest, for instance, can be 
obtained by using the large baselines available with the 
interferometer experiments. This combined with observa- 
tional data on foregrounds at frequencies slightly higher 
than 1421 MHz should provide us with an improved fore- 
ground model prior to the data analysis. However, because 
of the huge foreground amplitude compared to the 21cm 
signal, even very small variations to the assumed model 
might be confused with the signal and bias the estimation. 
In particular, as we have seen above, a good knowledge of 
the foreground frequency coherence structure in the range 
of interest will be crucial to the analysis. Therefore, we 
will now allow for free parameters in our foreground model 
and consider the case in which these parameters, together 
with the 21cm power spectra are simultaneously estimated 
from the data. Our foreground model is as given in equa- 
tions 14 with 

-2) =(^l + ^^^j . (36) 

The parameters to estimate are p = {A, /?, a, ^, log(n)} 
for each of the four foregrounds we are considering (point 
sources, Galactic synchrotron, extragalactic free-free and 
Galactic free-free emission) and Pa = Cf^{va, ^a)- We use 
49 bins in I space at each frequency (from I = 200 to 
5000 with a; = 100) and 40 bins of IMHz from 120MHz 
to 160MHz (z = 7.9 to z = 10.8) giving a total of 1980 
parameters (20 parameters for the foregrounds plus 1960 
for Cf^{va,Va))- Each bin in / space has Ni = 



independent aim measurements so that there will be a total 
of 642605 input data points in the form of aim values. 
We show in Figure 13 the estimated error on the power 
spectrum for the frequency bin of 140MHz from our full 
calculation. Using these 40 frequency bins, the foreground 
removal is quite efficient with the errors only showing a 
small foreground contamination. Note however that these 
errors will increase as we move toward the limits of the 
frequency range, as shown earlier in Figure 9. 

6. CONSTRAINING THE REIONIZATION MODEL 

Given a model of the perturbations in the ionization 
fraction we can ask how well one can constrain the corre- 
sponding free parameters. Typically, the ionization frac- 
tion power spectrum can be divided into two regions. First, 
on large enough scales, we expect the fraction of ionized 
material to follow the number of halos above a certain 
mass. The power spectrum will then be proportional to 
the matter density, with a bias factor that is expected to be 
larger than one. Second, as the scales approach the ionized 
bubble size, the power should be smoothed out since there 
are no fluctuations from inside the bubbles (this small scale 
cutoff might be quite extended if there is a wide distribu- 
tion of bubble sizes). The transition from one region to 
the other is model dependent. We are probing redshifts 
from z 8 to z 11 up to / ~ 1000 corresponding to 
scales X > 9Mpc which we expect to be larger than the 
typical bubble size in our model. We therefore model the 
ionization fraction power spectrum as 

Ps^sAk,z)^b^z)Pss{k,z) (37) 

and choose our free parameter to be bx = Xeb{z). For 
large bias (6(z)), the main contribution to equation 11 is 
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Fig. 12. — The recovered error at I = 1000 and u = 140 MHz, as- 
suming point sources contamination and using the correlation func- 
tions from figure 11 (blue solid line). Dashed line shows the signal. 



Fig. 13. — Expected error on the power spectrum at 140MHz, 
when including point sources, Galactic synchrotron, extragalactic 
free- free and Galactic free-free (blue solid line). Dashed line gives 
the case with no foregrounds. Dot-dashed line shows the signal. 
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then b1{z)Pss{k, z). We have therefore a total of 60 pa- 
rameters in our model (one bx for each of the 40 frequency 
bins plus 20 foreground parameters). Figure 14 shows the 
expected error from our full calculation. For comparison 
we also plot the expected bx from two other reionization 
models with the same optical depth (see Figures 3 and 
4). The recovered error is only at a few percent level so 
that, by measuring this parameter, it should be possible 
to distinguish between different reionization histories. 

7. SUMMARY 

The measurement of the 21cm radiation from neutral 
hydrogen at high redshifts will provide unprecedent infor- 
mation about the epoch of reionization. Unfortunately 
this signal has a huge foreground contamination, several 
orders of magnitude above the 21cm fluctuations (Oh & 
Mack 2003), making its detection a real challenge. In 
this paper we have presented a detailed study of the fore- 
grounds and the manner in which they will affect the mea- 
surement of the 21cm perturbations. We consider four 
types of foregrounds: diffuse Galactic synchrotron emis- 
sion at low radio frequencies. Galactic free-free, extra- 
galactic free-free emission and extragalactic point sources. 
We analyzed both the case when the foreground power 
spectra are known and the case when their parameters 
must be determined from the 21cm data itself. 

Our results are in general encouraging and show that 
we can achieve a high level of foreground cleaning even 
when using quite conservative choices for foreground mod- 
els both in terms of foreground fluctuation amplitude and 
frequency coherence. The removal is aided by the large fre- 
quency cross-correlation for foreground signals, although 
their frequency dependence and spatial power spectra are 
also useful in distinguishing them from the 21 cm signal. 
The 21 cm signal correlations, between channels, also play 
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Fig. 14. — Expected error on the product of the bias times the 
ionization fraction. Models are the same as in Figures 3 and 4. 



an important role and basically set a minimum bandwidth 
to the frequency cleaning process in which foreground re- 
moval is still effective. We do not expect a binning finer 
than IMHz to greatly improve the reconstruction. Also, as 
expected, foreground removal improves as we increase the 
number of frequency bins with fixed bin size. Using the 
full information from 40 bins of IMHz each, it is possible 
to clean out the foreground signal almost completely from 
the background 21 cm signal. 

The most crucial feature in the cleaning process is the 
foreground correlation structure. Concentrating on the 
point source contamination (which should be the most 
damaging foreground due to its smaller coherence), we see 
that the Gaussian correlation function provides the best 
results. In this case, limits on the spectral index variation 
are well above expected values from point source mod- 
els. Even in the worst case scenario, when the function 
is nearly linear in log (vi / 1^2) , foreground removal is still 
effective as long as the spectral index variations Aa < 1, 
which is consistent with observations of point sources (Co- 
hen et al. 2004). Our full statistical analysis allowed for 
variations in this correlation structure as well as the spec- 
tral index and shape of the angular power spectrum. Using 
a total of 40 bins from 120MHz to 160MHz and taking into 
account all the foregrounds, we calculated the expected er- 
rors for a total of 1975 parameters. The results show that 
foreground contamination will be minimal, allowing for a 
clear measurement of 21cm power spectra. 

We also investigated how well one can constrain the 
reionization model, by making full use of the three dimen- 
sional nature of the data (both in frequency and angular 
dependence) . Although a detailed parameterization of the 
reionization model incorporating all possible variations to 
the reionization history is difficult, we expect certain fea- 
tures to be robust. On scales larger than the patch size, for 
instance, the 21cm power spectra should follow the shape 
of the projected dark matter power spectrum. Our analy- 
sis has shown that the large-scale bias of the neutral Hy- 
drogen gas distribution, with respect to the density field, 
can be determined with a high precision. This allows for 
an accurate distinction between certain reionization mod- 
els, even if they have the same optical depth. 

The technical requirements to measure the 21cm signal, 
although challenging, should be within reach by proposed 
experiments such as SKA, LOFAR and PAST. Instruments 
such as PAST can provide additional information related 
to foregrounds that can be used to optimize SKA/LOFAR 
surveys. The big question addressed in this paper was 
to quantitatively establish if foregrounds, when taken as 
a whole, would contaminate 21 cm measurements beyond 
repair. Our results show that the future for these exper- 
iments is quite promising: even fully taking into account 
the expected foregrounds, the 21 cm fluctuations can be 
measured well enough to be a powerful probe of reioniza- 
tion. 
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